CAPTAIN2 workflow

EDGE2 continental 0.1 budget

library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(smoothr)
library(raster)
library(readr)
library(tidyr)

# Read the CAPTAIN2 EDGE2 RDS file
CAPTAIN2_EDGE2_data <- readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))

# Analyze non-zero cells in the prioritization output
cat("Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:\n")

Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:

total_cells <- nrow(CAPTAIN2_EDGE2_data)
nonzero_cells <- sum(CAPTAIN2_EDGE2_data$Priority > 0, na.rm = TRUE)
zero_cells <- sum(CAPTAIN2_EDGE2_data$Priority == 0, na.rm = TRUE)
na_cells <- sum(is.na(CAPTAIN2_EDGE2_data$Priority))

cat("Total cells in grid:", total_cells, "\n")

Total cells in grid: 232560 

cat("Cells with non-zero priority:", nonzero_cells, " (", round(nonzero_cells/total_cells*100, 2), "%)\n", sep="")

Cells with non-zero priority:5618 (2.42%)

cat("Cells with zero priority:", zero_cells, " (", round(zero_cells/total_cells*100, 2), "%)\n", sep="")

Cells with zero priority:226942 (97.58%)

cat("Cells with NA priority:", na_cells, " (", round(na_cells/total_cells*100, 2), "%)\n", sep="")

Cells with NA priority:0 (0%)

# Summary statistics of priority values
priority_summary <- summary(CAPTAIN2_EDGE2_data$Priority)
cat("\nSummary statistics of priority values:\n")

Summary statistics of priority values:

print(priority_summary)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.00000 0.00000 0.02331 0.00000 1.00000 

# Distribution of non-zero priority values
nonzero_priority <- CAPTAIN2_EDGE2_data$Priority[CAPTAIN2_EDGE2_data$Priority > 0]
cat("\nDistribution of non-zero priority values:\n")

Distribution of non-zero priority values:

priority_quantiles <- quantile(nonzero_priority, probs = seq(0, 1, 0.1), na.rm = TRUE)
print(priority_quantiles)

  0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
0.02 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 

# Load one of your input raster files to extract the correct grid structure
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 data based on PUID
CAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_EDGE2_sf <- st_as_sf(
  CAPTAIN2_EDGE2_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_EDGE2 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_EDGE2 <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_EDGE2 <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# Create the plot
CAPTAIN2_EDGE2_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the plot
print(CAPTAIN2_EDGE2_plot)

# Save the plot
ggsave(
  filename = here::here("outputs", "CAPTAIN2_2ndrun_EDGE2_priorities_01_2.png"),
  plot = CAPTAIN2_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

FUSE continental 0.1 budget

library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(smoothr)
library(raster)

# Read the CAPTAIN2 FUSE RDS file
CAPTAIN2_FUSE_data <- readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))

# Load one of your input raster files to extract the correct grid structure
# Use the same raster file for consistency
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 FUSE data based on PUID
CAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_FUSE_sf <- st_as_sf(
  CAPTAIN2_FUSE_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_FUSE <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_FUSE <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_FUSE <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# Create the plot
CAPTAIN2_FUSE_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE

# Display the plot
print(CAPTAIN2_FUSE_plot)

# Save the plot
ggsave(
  filename = here::here("outputs", "CAPTAIN2_FUSE_priorities_01_2.png"),
  plot = CAPTAIN2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

IUCN continental 0.1 budget

library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(smoothr)
library(raster)

# Read the CAPTAIN2 IUCN RDS file
CAPTAIN2_IUCN_data <- readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))

# Load one of your input raster files to extract the correct grid structure
# Use the same raster file for consistency
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Get the dimensions of the raster
nrows <- nrow(r)
ncols <- ncol(r)

# Confirm dimensions match expected values
if (nrows != 323 || ncols != 720) {
  warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")
}

# Create a grid of coordinates for each cell
# This gives us the center coordinates of each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Now join with the CAPTAIN2 IUCN data based on PUID
CAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%
  left_join(coords, by = "PUID")

# Check if the join worked correctly
if (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) > 0) {
  warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")
}

# Filter to keep only cells with non-zero priority for faster plotting
CAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%
  filter(Priority > 0) %>%
  filter(!is.na(Longitude), !is.na(Latitude))  # Remove any rows with missing coords

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform the dataset to sf object and project
CAPTAIN2_IUCN_sf <- st_as_sf(
  CAPTAIN2_IUCN_data_nonzero, 
  coords = c("Longitude", "Latitude"), 
  crs = crs(r, asText = TRUE)  # Use the raster's CRS
) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected_CAPTAIN2_IUCN <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border_CAPTAIN2_IUCN <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme_CAPTAIN2_IUCN <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# Create the plot
CAPTAIN2_IUCN_plot <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Global Conservation Priorities",
       subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN

# Display the plot
print(CAPTAIN2_IUCN_plot)

# Save the plot
ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_priorities_01_2.png"),
  plot = CAPTAIN2_IUCN_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

Difference maps

library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(smoothr)
library(raster)

# Read all three index RDS files
CAPTAIN2_IUCN_data <- readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))
CAPTAIN2_EDGE2_data <- readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))
CAPTAIN2_FUSE_data <- readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))

# Load one of your input raster files to extract the correct grid structure
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")

# Check if the file exists
if (!file.exists(raster_file)) {
  stop("Raster file not found. Please provide a valid path to one of your input raster files.")
}

# Load the raster
r <- raster(raster_file)

# Create a grid of coordinates for each cell
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")

# Add cell IDs (PUID) to the coordinates
coords$PUID <- 1:nrow(coords)

# Join all three datasets with coordinates
IUCN_with_coords <- CAPTAIN2_IUCN_data %>%
  dplyr::select(PUID, IUCN = Priority) %>%
  left_join(coords, by = "PUID")

EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>%
  dplyr::select(PUID, EDGE2 = Priority) %>%
  left_join(coords, by = "PUID") 

FUSE_with_coords <- CAPTAIN2_FUSE_data %>%
  dplyr::select(PUID, FUSE = Priority) %>%
  left_join(coords, by = "PUID")

# Combine all datasets
all_indices <- coords %>%
  left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by = "PUID")

# Calculate differences
all_indices <- all_indices %>%
  mutate(
    # Replace NA with 0 for calculation purposes
    IUCN = ifelse(is.na(IUCN), 0, IUCN),
    EDGE2 = ifelse(is.na(EDGE2), 0, EDGE2),
    FUSE = ifelse(is.na(FUSE), 0, FUSE),
    
    # Calculate differences
    IUCN_minus_FUSE = IUCN - FUSE,
    IUCN_minus_EDGE2 = IUCN - EDGE2,
    EDGE2_minus_FUSE = EDGE2 - FUSE
  )

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# Filter to non-zero differences for each comparison to reduce plot size
# IUCN - FUSE
IUCN_FUSE_diff <- all_indices %>%
  filter(IUCN_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# IUCN - EDGE2
IUCN_EDGE2_diff <- all_indices %>%
  filter(IUCN_minus_EDGE2 != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# EDGE2 - FUSE
EDGE2_FUSE_diff <- all_indices %>%
  filter(EDGE2_minus_FUSE != 0) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(r, asText = TRUE)) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create a diverging color palette for difference maps
# Blue for negative (first index lower), white for zero, red for positive (first index higher)
diff_colors <- c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", 
                "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")

# 1. IUCN - FUSE Difference Map
IUCN_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority\n(IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority\n(IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_plot <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority\n(EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Difference in Conservation Priorities",
       subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Display all plots
print(IUCN_FUSE_plot)

print(IUCN_EDGE2_plot)

print(EDGE2_FUSE_plot)

# Save all plots
ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_FUSE_difference_2.png"),
  plot = IUCN_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_IUCN_minus_EDGE2_difference_2.png"),
  plot = IUCN_EDGE2_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

ggsave(
  filename = here::here("outputs", "CAPTAIN2_EDGE2_minus_FUSE_difference_2.png"),
  plot = EDGE2_FUSE_plot,
  width = 10,
  height = 6,
  dpi = 300,
  bg = "white"
)

# Optionally, create a panel with all three difference maps
library(patchwork)

# Combine all plots
all_diffs_plot <- IUCN_FUSE_plot / IUCN_EDGE2_plot / EDGE2_FUSE_plot +
  plot_annotation(
    title = "Differences Between Conservation Priority Indices",
    subtitle = "Budget: 0.1, Replicates: 100",
    theme = theme(plot.title = element_text(hjust = 0.5),
                  plot.subtitle = element_text(hjust = 0.5))
  )

#all_diffs_plot

# Save the combined plot
#ggsave(
#  filename = here::here("outputs", "CAPTAIN2_all_differences.png"),
#  plot = all_diffs_plot,
#  width = 10,
#  height = 15,
#  dpi = 300,
#  bg = "white"
#)

Species level priorities

# Load required packages
library(here)
library(readr)
library(dplyr)
library(tidyr)
library(ggplot2)

# Read the protected range fractions RDS file
protected_fractions <- readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_2ndrun.rds"))

# Read the continental shark conservation metrics CSV file
shark_metrics <- read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))

# Rename column in shark_metrics to match better
shark_metrics <- shark_metrics %>%
  rename(Species = `Species name`)

# Order shark_metrics alphabetically by Species name
shark_metrics <- shark_metrics %>%
  arrange(Species)

# Add an original order ID to protected_fractions to maintain its original order
protected_fractions$original_order <- 1:nrow(protected_fractions)

# Add row number as IDs to both datasets
protected_fractions$protected_ID <- 1:nrow(protected_fractions)
shark_metrics$species_ID <- 1:nrow(shark_metrics)

# Check if the datasets have the same number of rows
if(nrow(protected_fractions) == nrow(shark_metrics)) {
  # Create index mapping - this maintains the original protection data ordering
  # while allowing us to associate with alphabetically ordered species names
  indices <- data.frame(
    protected_ID = 1:nrow(protected_fractions),
    species_ID = 1:nrow(shark_metrics)
  )
  
  # Join protected_fractions with indices
  protected_with_indices <- protected_fractions %>%
    left_join(indices, by = "protected_ID")
  
  # Join shark_metrics with indices
  shark_with_indices <- shark_metrics %>%
    left_join(indices, by = "species_ID")
  
  # Now join the datasets, matching on species_ID and protected_ID
  combined_data <- protected_with_indices %>%
    inner_join(
      shark_with_indices,
      by = c("species_ID", "protected_ID"),
      suffix = c("_captain", "_original")
    ) %>%
    # Sort by the original order of protected_fractions 
    arrange(original_order)
  
  cat("Successfully joined datasets with", nrow(combined_data), "species\n")
  cat("First few species in combined dataset:\n")
  print(head(combined_data[, c("Species_captain", "Species_original")]))
  
  # Define IUCN categories and order - using only the first 5 categories
  iucn_labels <- c(
    "1" = "LC", 
    "2" = "NT", 
    "3" = "VU", 
    "4" = "EN", 
    "5" = "CR"
  )
  
  iucn_order <- c("LC", "NT", "VU", "EN", "CR")
  
  # Define colors for IUCN categories
  iucn_colors <- c(
    "LC" = "#50C878",     # Green
    "NT" = "#FFFF00",     # Yellow
    "VU" = "#FFA500",     # Orange
    "EN" = "#FF8C00",     # Dark Orange
    "CR" = "#FF0000"      # Red
  )
  
  # 1. IUCN Boxplot
  iucn_boxplot <- combined_data %>%
    mutate(
      IUCN_status = factor(IUCN_original, levels = 1:5, labels = iucn_order),
      protection_percentage = IUCN_captain * 100
    ) %>%
    ggplot(aes(x = IUCN_status, y = protection_percentage)) +
  #  geom_violin(aes(fill = IUCN_status, color = IUCN_status), 
  #              trim = FALSE, 
  #              alpha = 0.5) +
    geom_jitter(width = 0.1, 
                size = 0.6, 
                alpha = 0.5, 
                color = "darkgray") +
    geom_boxplot(width = 0.1, 
                 fill = "white", 
                 color = "black", 
                 outlier.shape = NA, 
                 alpha = 0.8) +
    labs(title = "IUCN Priority Index: Range Protection by IUCN Status",
         x = "IUCN Red List threat status", 
         y = "Range protected (%)") +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      legend.position = "none",
      panel.grid.major.x = element_blank()
    ) +
    scale_fill_manual(values = iucn_colors) +
    scale_color_manual(values = iucn_colors) +
    scale_y_continuous(limits = c(0, 100),
                      breaks = seq(0, 100, 25))
  
  # 2. FUSE Boxplot
  fuse_boxplot <- combined_data %>%
    mutate(
      FUSE_category = factor(FUSE_original),
      protection_percentage = FUSE_captain * 100
    ) %>%
    ggplot(aes(x = FUSE_category, y = protection_percentage)) +
  #  geom_violin(aes(fill = FUSE_category, color = FUSE_category), 
  #              trim = FALSE, 
  #              alpha = 0.5) +
    geom_jitter(width = 0.1, 
                size = 0.6, 
                alpha = 0.5, 
                color = "darkgray") +
    geom_boxplot(width = 0.1, 
                 fill = "white", 
                 color = "black", 
                 outlier.shape = NA, 
                 alpha = 0.8) +
    labs(title = "FUSE Priority Index: Range Protection by FUSE Score",
         x = "FUSE Score", 
         y = "Range protected (%)") +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      legend.position = "none",
      panel.grid.major.x = element_blank()
    ) +
    scale_y_continuous(limits = c(0, 100),
                      breaks = seq(0, 100, 25))
  
  # 3. EDGE2 Boxplot
  edge2_boxplot <- combined_data %>%
    mutate(
      EDGE2_category = factor(EDGE2_original),
      protection_percentage = EDGE2_captain * 100
    ) %>%
    ggplot(aes(x = EDGE2_category, y = protection_percentage)) +
  #  geom_violin(aes(fill = EDGE2_category, color = EDGE2_category), 
  #              trim = FALSE, 
  #              alpha = 0.5) +
    geom_jitter(width = 0.1, 
                size = 0.6, 
                alpha = 0.5, 
                color = "darkgray") +
    geom_boxplot(width = 0.1, 
                 fill = "white", 
                 color = "black", 
                 outlier.shape = NA, 
                 alpha = 0.8) +
    labs(title = "EDGE2 Priority Index: Range Protection by EDGE2 Score",
         x = "EDGE2 Score", 
         y = "Range protected (%)") +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      legend.position = "none",
      panel.grid.major.x = element_blank()
    ) +
    scale_y_continuous(limits = c(0, 100),
                      breaks = seq(0, 100, 25))
  
  # Display individual plots
  print(iucn_boxplot)
  print(fuse_boxplot)
  print(edge2_boxplot)
  
  # Save individual plots
  ggsave(here::here("outputs", "iucn_protection_boxplot.png"), 
         iucn_boxplot, width = 10, height = 6, dpi = 300)
  ggsave(here::here("outputs", "fuse_protection_boxplot.png"), 
         fuse_boxplot, width = 10, height = 6, dpi = 300)
  ggsave(here::here("outputs", "edge2_protection_boxplot.png"), 
         edge2_boxplot, width = 10, height = 6, dpi = 300)
  
  # Create a nicely formatted table showing the ordered species
  species_table <- combined_data %>%
    dplyr::select(Species_captain, Species_original, IUCN_captain, FUSE_captain, EDGE2_captain, 
           IUCN_original, FUSE_original, EDGE2_original) %>%
    arrange(Species_original) %>%
    head(20)  # Just show the first 20 for display
  
  # Print species table
  cat("\nFirst 20 species (alphabetically by original species name):\n")
  print(species_table)
  
  # Save the full combined data
  write.csv(combined_data, here::here("outputs", "combined_species_data.csv"), row.names = FALSE)
  
} else {
  cat("ERROR: Datasets have different number of rows.\n")
  cat("Protected fractions:", nrow(protected_fractions), "rows\n")
  cat("Shark metrics:", nrow(shark_metrics), "rows\n")
}

Successfully joined datasets with 1000 species
First few species in combined dataset:
  Species_captain           Species_original
1       Species_1   Acroteriobatus_annulatus
2       Species_2     Acroteriobatus_blochii
3       Species_3 Acroteriobatus_leucospilus
4       Species_4   Acroteriobatus_ocellatus
5       Species_5     Acroteriobatus_salalah
6       Species_6  Acroteriobatus_variegatus

First 20 species (alphabetically by original species name):
   Species_captain           Species_original IUCN_captain FUSE_captain
1        Species_1   Acroteriobatus_annulatus   0.88645161 0.5374193548
2        Species_2     Acroteriobatus_blochii   0.18000000 0.5700000000
3        Species_3 Acroteriobatus_leucospilus   0.87462500 0.6245000000
4        Species_4   Acroteriobatus_ocellatus   1.00000000 0.6711111111
5        Species_5     Acroteriobatus_salalah   0.54017699 0.1313274336
6        Species_6  Acroteriobatus_variegatus   0.84383562 0.5056164384
7        Species_7 Acroteriobatus_zanzibarens   0.76666667 0.0600000000
8        Species_8             Aculeola_nigra   0.37897436 0.0005128205
9        Species_9        Aetobatus_flagellum   0.83666667 0.4534920635
10      Species_10         Aetobatus_narinari   0.32685917 0.2346902159
11      Species_11     Aetomylaeus_asperrimus   1.00000000 0.2380952381
12      Species_12        Aetomylaeus_bovinus   0.34552910 0.7692063492
13      Species_13      Aetomylaeus_maculatus   0.99058366 0.9174319066
14      Species_14       Aetomylaeus_nichofii   0.44169283 0.4223268698
15      Species_15    Aetomylaeus_vespertilio   0.76333492 0.4777227251
16      Species_16          Alopias_pelagicus   0.11931061 0.0890257673
17      Species_17      Alopias_superciliosus   0.14376233 0.1325789344
18      Species_18           Alopias_vulpinus   0.12449208 0.1288068917
19      Species_19    Amblyraja_doellojuradoi   0.31948148 0.0146296296
20      Species_20       Amblyraja_hyperborea   0.04166804 0.1077711340
   EDGE2_captain IUCN_original FUSE_original EDGE2_original
1     0.00000000             3             1              1
2     0.00000000             1             1              1
3     0.02775000             4             1              1
4     0.07851852             5             1              2
5     0.02761062             2             1              1
6     0.44671233             5             2              2
7     0.00000000             1             1              1
8     1.00000000             2             1              1
9     0.24047619             4             3              3
10    0.24811700             4             2              3
11    1.00000000             1             1              1
12    0.77567460             1             4              1
13    0.90428016             4             4              2
14    0.29180055             3             1              1
15    0.57068128             4             2              2
16    0.07715753             4             2              2
17    0.13980821             3             1              1
18    0.12752863             3             1              1
19    0.47288889             1             1              1
20    0.09284124             1             1              1

Manuscript maps

Individual maps

CAPTAIN2_EDGE2_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2
CAPTAIN2_EDGE2_msmap

CAPTAIN2_FUSE_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_FUSE, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_FUSE, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_FUSE
CAPTAIN2_FUSE_msmap

CAPTAIN2_IUCN_msmap <- ggplot() +
  geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_CAPTAIN2_IUCN, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_IUCN, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Global Conservation Priorities",
       #subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_IUCN
CAPTAIN2_IUCN_msmap

# Load required libraries
library(patchwork)
library(ggplot2)

# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
CAPTAIN2_IUCN_msmap_labeled <- add_panel_labels(CAPTAIN2_IUCN_msmap, "A")
CAPTAIN2_FUSE_msmap_labeled <- add_panel_labels(CAPTAIN2_FUSE_msmap, "B")
CAPTAIN2_EDGE2_msmap_labeled <- add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid1 <- CAPTAIN2_IUCN_msmap_labeled /
         CAPTAIN2_FUSE_msmap_labeled /
         CAPTAIN2_EDGE2_msmap_labeled

# Display each grid separately
grid1

# Save the grids if needed
ggsave(here::here("grid1_maps_continental_2ndrun.png"), grid1, width = 8, height = 15, dpi = 300)

Difference maps

# 1. IUCN - FUSE Difference Map
IUCN_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority\n(IUCN - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# 2. IUCN - EDGE2 Difference Map
IUCN_EDGE2_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority\n(IUCN - EDGE2)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "IUCN Index minus EDGE2 Index",
       x = NULL, y = NULL) +
  my_theme

# 3. EDGE2 - FUSE Difference Map
EDGE2_FUSE_msmap <- ggplot() +
  geom_sf(data = globe_border, fill = "#F8F8F8", color = NA) +
  geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size = 0.5) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "darkgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = diff_colors,
    limits = c(-1, 1),
    breaks = seq(-1, 1, by = 0.25),
    labels = as.character(seq(-1, 1, by = 0.25)),
    name = "Difference in Priority\n(EDGE2 - FUSE)",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(#title = "Difference in Conservation Priorities",
       #subtitle = "EDGE2 Index minus FUSE Index",
       x = NULL, y = NULL) +
  my_theme

# Load required libraries
library(patchwork)
library(ggplot2)

# Create a function to add labels (A), (B), etc.
add_panel_labels <- function(plot, label) {
  plot + 
    theme(
      plot.title = element_text(face = "bold", hjust = 0, size = 12)
    ) +
    labs(title = paste0("(", label, ")"))
}

# Add labels to each plot
# First grid
IUCN_FUSE_msmap_labeled <- add_panel_labels(IUCN_FUSE_msmap, "A")
IUCN_EDGE2_msmap_labeled <- add_panel_labels(IUCN_EDGE2_msmap, "B")
EDGE2_FUSE_msmap_labeled <- add_panel_labels(EDGE2_FUSE_msmap, "C")

# Combine plots into two separate grids, each with 3 rows and 1 column
grid2 <- IUCN_FUSE_msmap_labeled /
         IUCN_EDGE2_msmap_labeled /
         EDGE2_FUSE_msmap_labeled

# Display each grid separately
grid2

# Save the grids if needed
ggsave(here::here("grid2_maps_continental_2ndrun.png"), grid2, width = 8, height = 15, dpi = 300)

Congruence maps

# Method 1: If you have separate dataframes for each index
# Assuming your data has columns: PUID, Priority, and geometry

# First, identify high priority areas (>0.9) for each index
high_priority_EDGE2 <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$Priority > 0.9, ]
high_priority_FUSE <- CAPTAIN2_FUSE_sf[CAPTAIN2_FUSE_sf$Priority > 0.9, ]
high_priority_IUCN <- CAPTAIN2_IUCN_sf[CAPTAIN2_IUCN_sf$Priority > 0.9, ]

# Find congruent areas (present in all three)
# Method using PUID (assuming you have PUID column)
congruent_PUIDs <- intersect(intersect(high_priority_EDGE2$PUID, 
                                      high_priority_FUSE$PUID), 
                            high_priority_IUCN$PUID)

# Extract congruent areas
congruent_areas <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$PUID %in% congruent_PUIDs, ]

# Print summary
cat("High priority areas (>0.9):\n")

High priority areas (>0.9):

cat("EDGE2:", nrow(high_priority_EDGE2), "areas\n")

EDGE2: 5328 areas

cat("FUSE:", nrow(high_priority_FUSE), "areas\n") 

FUSE: 4396 areas

cat("IUCN:", nrow(high_priority_IUCN), "areas\n")

IUCN: 4734 areas

cat("Congruent areas:", nrow(congruent_areas), "areas\n")

Congruent areas: 1459 areas

cat("Percentage of overlap:", round(nrow(congruent_areas)/min(nrow(high_priority_EDGE2), 
                                                            nrow(high_priority_FUSE), 
                                                            nrow(high_priority_IUCN))*100, 2), "%\n")

Percentage of overlap: 33.19 %

# Method 2: If you need to merge dataframes first
# Create a combined dataframe with all three priority scores
# First, extract just the data (without geometry) from the other sf objects
FUSE_data <- st_drop_geometry(CAPTAIN2_FUSE_sf[, c("PUID", "Priority")])
IUCN_data <- st_drop_geometry(CAPTAIN2_IUCN_sf[, c("PUID", "Priority")])

# Merge with the EDGE2 sf object (keeping geometry)
combined_priorities <- merge(CAPTAIN2_EDGE2_sf, FUSE_data, 
                            by = "PUID", suffixes = c("_EDGE2", "_FUSE"))
combined_priorities <- merge(combined_priorities, IUCN_data, 
                            by = "PUID")
names(combined_priorities)[names(combined_priorities) == "Priority"] <- "Priority_IUCN"

# Identify congruent high priority areas
congruent_areas_v2 <- combined_priorities[combined_priorities$Priority_EDGE2 > 0.9 & 
                                         combined_priorities$Priority_FUSE > 0.9 & 
                                         combined_priorities$Priority_IUCN > 0.9, ]

# Create a map showing only congruent areas
congruent_map <- ggplot() +
  geom_sf(data = congruent_areas, aes(fill = "Congruent High Priority"), 
          color = "red", size = 0.8, alpha = 0.8) +
  geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill = NA, color = "black", size = 0.5) +
  scale_fill_manual(values = c("Congruent High Priority" = "red"),
                    name = "Priority Areas") +
  labs(title = "Congruent High Priority Areas",
       subtitle = "Areas with Priority > 0.9 in EDGE2, FUSE, and IUCN indices",
       x = NULL, y = NULL) +
  my_theme_CAPTAIN2_EDGE2

# Display the map
print(congruent_map)

# Create detailed congruence categories
# Define logical conditions for each index
EDGE2_high <- combined_priorities$Priority_EDGE2 > 0.9
FUSE_high <- combined_priorities$Priority_FUSE > 0.9
IUCN_high <- combined_priorities$Priority_IUCN > 0.9

# Create detailed congruence categories
combined_priorities$congruence_category <- case_when(
  EDGE2_high & FUSE_high & IUCN_high ~ "All three indices",
  EDGE2_high & FUSE_high & !IUCN_high ~ "EDGE2 + FUSE",
  EDGE2_high & !FUSE_high & IUCN_high ~ "EDGE2 + IUCN", 
  !EDGE2_high & FUSE_high & IUCN_high ~ "FUSE + IUCN",
  EDGE2_high & !FUSE_high & !IUCN_high ~ "EDGE2 only",
  !EDGE2_high & FUSE_high & !IUCN_high ~ "FUSE only",
  !EDGE2_high & !FUSE_high & IUCN_high ~ "IUCN only",
  TRUE ~ "No high priority"
)

# Convert to factor with desired order
combined_priorities$congruence_category <- factor(combined_priorities$congruence_category,
                                                 levels = c("All three indices", 
                                                           "EDGE2 + FUSE", "EDGE2 + IUCN", "FUSE + IUCN",
                                                           "EDGE2 only", "FUSE only", "IUCN only"))

# Filter data to only include high priority areas
high_priority_data <- combined_priorities[combined_priorities$congruence_category %in% 
                                         c("All three indices", "EDGE2 + FUSE", "EDGE2 + IUCN", 
                                           "FUSE + IUCN", "EDGE2 only", "FUSE only", "IUCN only"), ]

# Debug: Check the data
cat("Data check:\n")

Data check:

cat("Total high priority areas:", nrow(high_priority_data), "\n")

Total high priority areas: 3139 

print(table(high_priority_data$congruence_category))

All three indices      EDGE2 + FUSE      EDGE2 + IUCN       FUSE + IUCN 
             1459               275              1130                75 
       EDGE2 only         FUSE only         IUCN only 
              154                17                29 

# Check if the factor levels are properly set
cat("\nFactor levels:\n")

Factor levels:

print(levels(high_priority_data$congruence_category))

[1] "All three indices" "EDGE2 + FUSE"      "EDGE2 + IUCN"     
[4] "FUSE + IUCN"       "EDGE2 only"        "FUSE only"        
[7] "IUCN only"        

# Check for any issues with the data
cat("\nData structure check:\n")

Data structure check:

print(str(high_priority_data$congruence_category))

 Factor w/ 7 levels "All three indices",..: 6 6 6 6 6 6 6 6 2 2 ...
NULL

# Test a simple plot first
cat("\nTesting simple plot...\n")

Testing simple plot...

test_plot <- ggplot() +
  geom_sf(data = high_priority_data, aes(fill = congruence_category)) +
  scale_fill_manual(values = c("All three indices" = "#8B0000",
                              "EDGE2 + FUSE" = "#FF4500",
                              "EDGE2 + IUCN" = "#FF6347",
                              "FUSE + IUCN" = "#FFA500",
                              "EDGE2 only" = "#4169E1",
                              "FUSE only" = "#32CD32",
                              "IUCN only" = "#9370DB")) +
  theme_void()

# Try to print the simple test plot
tryCatch({
  print(test_plot)
  cat("Simple plot works!\n")
}, error = function(e) {
  cat("Simple plot failed with error:", e$message, "\n")
})

Simple plot works!

# If simple plot fails, let's check the data more thoroughly
if (nrow(high_priority_data) == 0) {
  cat("No high priority data found! Checking original data...\n")
  cat("EDGE2 > 0.9:", sum(combined_priorities$Priority_EDGE2 > 0.9, na.rm = TRUE), "\n")
  cat("FUSE > 0.9:", sum(combined_priorities$Priority_FUSE > 0.9, na.rm = TRUE), "\n") 
  cat("IUCN > 0.9:", sum(combined_priorities$Priority_IUCN > 0.9, na.rm = TRUE), "\n")
}

# Extract coordinates for the congruence plot
if (nrow(high_priority_data) > 0) {
  coords <- st_coordinates(st_centroid(high_priority_data))
  plot_data <- data.frame(
    x = coords[,1],
    y = coords[,2], 
    category = high_priority_data$congruence_category
  )
  
  # Remove any rows with NA category
  plot_data <- plot_data[!is.na(plot_data$category), ]
  
  cat("Creating enhanced congruence plot...\n")
  cat("Plot data dimensions:", nrow(plot_data), "points\n")
  
  # Create the enhanced congruence plot
  congruence_map <- ggplot(plot_data, aes(x = x, y = y, color = category)) +
    geom_point(size = 1.2, alpha = 0.85, stroke = 0) +  # Larger points, no stroke to reduce overlap
    scale_color_manual(
      values = c("All three indices" = "#8B0000",      # Dark red
                "EDGE2 + FUSE" = "#FF8C00",            # Dark orange  
                "EDGE2 + IUCN" = "#DC143C",            # Crimson
                "FUSE + IUCN" = "#FFD700",             # Gold
                "EDGE2 only" = "#4169E1",              # Royal blue
                "FUSE only" = "#32CD32",               # Lime green
                "IUCN only" = "#9370DB"),              # Medium purple
      name = "High Priority\nCongruence",
      guide = guide_legend(
        override.aes = list(size = 4, alpha = 1),  # Larger, more opaque legend points
        title.position = "top",
        title.hjust = 0.5,
        ncol = 2  # Two columns for legend to save space
      )) +
    labs(
      title = "Global Conservation Priority Congruence Analysis",
      subtitle = "High priority areas (>0.9) showing agreement patterns across EDGE2, FUSE, and IUCN indices",
      x = "Longitude", 
      y = "Latitude",
      caption = paste("Total high priority areas:", nrow(plot_data))
    ) +
    theme_minimal() +
    theme(
      # Plot aesthetics
      plot.background = element_rect(fill = "white", color = NA),
      panel.background = element_rect(fill = "#f8f9fa", color = NA),
      panel.grid.major = element_line(color = "white", size = 0.5, linetype = "solid"),
      panel.grid.minor = element_line(color = "white", size = 0.25, linetype = "solid"),
      
      # Text styling
      plot.title = element_text(size = 16, hjust = 0.5, face = "bold", 
                               margin = margin(b = 10)),
      plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray30",
                                  margin = margin(b = 20)),
      plot.caption = element_text(size = 10, color = "gray50", hjust = 1),
      
      # Axis styling
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(size = 10, color = "gray30"),
      axis.ticks = element_line(color = "gray50", size = 0.5),
      
      # Legend styling
      legend.position = "bottom",
      legend.background = element_rect(fill = "white", color = "gray80", size = 0.5),
      legend.margin = margin(t = 15, r = 10, b = 10, l = 10),
      legend.title = element_text(size = 12, face = "bold"),
      legend.text = element_text(size = 10),
      legend.key.size = unit(0.8, "cm"),
      
      # Panel border
      panel.border = element_rect(color = "gray70", fill = NA, size = 0.5)
    ) +
    # Add coordinate system and aspect ratio
    coord_fixed(ratio = 1.3) +  # Adjust ratio for better map appearance
    # Add subtle borders around plot area
    scale_x_continuous(expand = c(0.02, 0.02)) +
    scale_y_continuous(expand = c(0.02, 0.02))
  
  # Print the enhanced plot
  print(congruence_map)
  cat("Enhanced congruence plot created successfully!\n")
  
  # Print summary by category
  cat("\nBreakdown by congruence category:\n")
  category_counts <- table(plot_data$category)
  category_percentages <- round(prop.table(category_counts) * 100, 1)
  
  for(i in 1:length(category_counts)) {
    cat(sprintf("%-20s: %4d points (%4.1f%%)\n", 
                names(category_counts)[i], 
                category_counts[i], 
                category_percentages[i]))
  }
  
} else {
  cat("No valid high priority data to plot\n")
}

Creating enhanced congruence plot...
Plot data dimensions: 3139 points

Enhanced congruence plot created successfully!

Breakdown by congruence category:
All three indices   : 1459 points (46.5%)
EDGE2 + FUSE        :  275 points ( 8.8%)
EDGE2 + IUCN        : 1130 points (36.0%)
FUSE + IUCN         :   75 points ( 2.4%)
EDGE2 only          :  154 points ( 4.9%)
FUSE only           :   17 points ( 0.5%)
IUCN only           :   29 points ( 0.9%)

# Print summary of congruence patterns
cat("\nCongruence Pattern Summary:\n")

Congruence Pattern Summary:

congruence_summary <- table(combined_priorities$congruence_category, useNA = "ifany")
print(congruence_summary)

All three indices      EDGE2 + FUSE      EDGE2 + IUCN       FUSE + IUCN 
             1459               275              1130                75 
       EDGE2 only         FUSE only         IUCN only              <NA> 
              154                17                29                 6 

# Calculate percentages
total_high_priority <- sum(congruence_summary[names(congruence_summary) != "No high priority"], na.rm = TRUE)
congruence_percentages <- round(congruence_summary / total_high_priority * 100, 2)
cat("\nPercentages of high priority areas:\n")

Percentages of high priority areas:

print(congruence_percentages[names(congruence_percentages) != "No high priority"])

All three indices      EDGE2 + FUSE      EDGE2 + IUCN       FUSE + IUCN 
            46.48              8.76             36.00              2.39 
       EDGE2 only         FUSE only         IUCN only              <NA> 
             4.91              0.54              0.92                   

Relationship with fishing effort : bar plots per country

library(here)
library(dplyr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(raster)
library(tidyr)
library(ggrepel)
library(ggsci)   # For nice color palettes

# Read all priority data
CAPTAIN2_EDGE2_data <- readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))
CAPTAIN2_FUSE_data <- readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))
CAPTAIN2_IUCN_data <- readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))

# Load fishing data
load(here::here("Data", "Raw", "Predicted_Fishing_Hours_05Deg.Rdata"))

# Load marine ecoregion shapefile
meow_ecos <- st_read(here("Data", "Shapefiles", "meow_ecos", "meow_ecos.shp"), quiet = TRUE)

# Load raster to get coordinate system
raster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")
r <- raster(raster_file)

# Create coordinates for all grid cells
coords <- as.data.frame(coordinates(r))
names(coords) <- c("Longitude", "Latitude")
coords$PUID <- 1:nrow(coords)

#---------------------- FISHING DATA PROCESSING ----------------------#

# Convert aggregated_data to an sf object
fishing_sf <- aggregated_data %>%
  filter(!is.na(lon_05deg), !is.na(lat_05deg), !is.na(predicted_fishing_hours)) %>%
  st_as_sf(coords = c("lon_05deg", "lat_05deg"), crs = 4326)

# Make sure CRS matches
st_crs(fishing_sf) <- st_crs(meow_ecos)

# Process fishing data for ECOREGIONS
fishing_with_ecoregion <- st_join(fishing_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))

ecoregion_fishing_stats <- fishing_with_ecoregion %>%
  st_drop_geometry() %>%
  group_by(ECOREGION, REALM) %>%
  summarize(
    mean_fishing_hours = mean(predicted_fishing_hours, na.rm = TRUE),
    median_fishing_hours = median(predicted_fishing_hours, na.rm = TRUE),
    fishing_q05 = quantile(predicted_fishing_hours, 0.05, na.rm = TRUE),
    fishing_q95 = quantile(predicted_fishing_hours, 0.95, na.rm = TRUE),
    fishing_cells = n(),
    .groups = 'drop'
  ) %>%
  filter(!is.na(ECOREGION)) %>%
  arrange(desc(mean_fishing_hours))

#---------------------- PRIORITY DATA PROCESSING ----------------------#

# Combine all priority data with coordinates
combined_data <- CAPTAIN2_EDGE2_data %>%
  rename(Priority_EDGE2 = Priority) %>%
  left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, Priority) %>% rename(Priority_FUSE = Priority), by = "PUID") %>%
  left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, Priority) %>% rename(Priority_IUCN = Priority), by = "PUID") %>%
  left_join(coords, by = "PUID")

# Convert to sf object for each index
iucn_sf <- combined_data %>%
  filter(!is.na(Longitude), !is.na(Latitude), !is.na(Priority_IUCN)) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)

fuse_sf <- combined_data %>%
  filter(!is.na(Longitude), !is.na(Latitude), !is.na(Priority_FUSE)) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)

edge2_sf <- combined_data %>%
  filter(!is.na(Longitude), !is.na(Latitude), !is.na(Priority_EDGE2)) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)

# Make sure CRS matches
st_crs(iucn_sf) <- st_crs(meow_ecos)
st_crs(fuse_sf) <- st_crs(meow_ecos)
st_crs(edge2_sf) <- st_crs(meow_ecos)

# Process priority data for ECOREGIONS - IUCN
iucn_with_ecoregion <- st_join(iucn_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))
ecoregion_iucn_stats <- iucn_with_ecoregion %>%
  st_drop_geometry() %>%
  group_by(ECOREGION, REALM) %>%
  summarize(
    mean_priority = mean(Priority_IUCN, na.rm = TRUE),
    median_priority = median(Priority_IUCN, na.rm = TRUE),
    priority_q05 = quantile(Priority_IUCN, 0.05, na.rm = TRUE),
    priority_q95 = quantile(Priority_IUCN, 0.95, na.rm = TRUE),
    priority_cells = n(),
    .groups = 'drop'
  ) %>%
  filter(!is.na(ECOREGION)) %>%
  arrange(desc(mean_priority))

# Process priority data for ECOREGIONS - FUSE
fuse_with_ecoregion <- st_join(fuse_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))
ecoregion_fuse_stats <- fuse_with_ecoregion %>%
  st_drop_geometry() %>%
  group_by(ECOREGION, REALM) %>%
  summarize(
    mean_priority = mean(Priority_FUSE, na.rm = TRUE),
    median_priority = median(Priority_FUSE, na.rm = TRUE),
    priority_q05 = quantile(Priority_FUSE, 0.05, na.rm = TRUE),
    priority_q95 = quantile(Priority_FUSE, 0.95, na.rm = TRUE),
    priority_cells = n(),
    .groups = 'drop'
  ) %>%
  filter(!is.na(ECOREGION)) %>%
  arrange(desc(mean_priority))

# Process priority data for ECOREGIONS - EDGE2
edge2_with_ecoregion <- st_join(edge2_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))
ecoregion_edge2_stats <- edge2_with_ecoregion %>%
  st_drop_geometry() %>%
  group_by(ECOREGION, REALM) %>%
  summarize(
    mean_priority = mean(Priority_EDGE2, na.rm = TRUE),
    median_priority = median(Priority_EDGE2, na.rm = TRUE),
    priority_q05 = quantile(Priority_EDGE2, 0.05, na.rm = TRUE),
    priority_q95 = quantile(Priority_EDGE2, 0.95, na.rm = TRUE),
    priority_cells = n(),
    .groups = 'drop'
  ) %>%
  filter(!is.na(ECOREGION)) %>%
  arrange(desc(mean_priority))

#---------------------- MERGE DATASETS FOR SCATTERPLOTS ----------------------#

# Merge ECOREGION data for each index
ecoregion_combined_iucn <- inner_join(
  ecoregion_fishing_stats, 
  ecoregion_iucn_stats, 
  by = c("ECOREGION", "REALM")) %>%
  mutate(cell_ratio = fishing_cells / priority_cells)

ecoregion_combined_fuse <- inner_join(
  ecoregion_fishing_stats, 
  ecoregion_fuse_stats, 
  by = c("ECOREGION", "REALM")) %>%
  mutate(cell_ratio = fishing_cells / priority_cells)

ecoregion_combined_edge2 <- inner_join(
  ecoregion_fishing_stats, 
  ecoregion_edge2_stats, 
  by = c("ECOREGION", "REALM")) %>%
  mutate(cell_ratio = fishing_cells / priority_cells)

#---------------------- CREATE ECOREGION SCATTERPLOTS ----------------------#

# IUCN - Get top 5 fishing and top 5 priority ecoregions
top_fishing_iucn <- ecoregion_combined_iucn %>%
  arrange(desc(mean_fishing_hours)) %>%
  head(5) %>%
  pull(ECOREGION)

top_priority_iucn <- ecoregion_combined_iucn %>%
  arrange(desc(mean_priority)) %>%
  head(5) %>%
  pull(ECOREGION)

ecoregions_to_label_iucn <- unique(c(top_fishing_iucn, top_priority_iucn))
label_data_iucn <- ecoregion_combined_iucn %>%
  filter(ECOREGION %in% ecoregions_to_label_iucn)

# Create IUCN plot
iucn_ecoregion_plot <- ggplot(ecoregion_combined_iucn,
                             aes(x = mean_fishing_hours, y = mean_priority)) +
  geom_point(aes(size = fishing_cells, color = REALM), alpha = 0.4) +
  geom_label_repel(
    data = label_data_iucn,
    aes(label = ECOREGION),
    size = 3,
    max.overlaps = 10,
    box.padding = 0.5
  ) +
  scale_x_log10() +
  scale_size_continuous(name = "Number of Cells", range = c(1, 8)) +
  scale_color_manual(values = c("red", "blue", "darkgreen", "purple", "orange", "brown", 
                               "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +
  labs(title = "Relationship Between Fishing Pressure and Conservation Priority",
       subtitle = "By Ecoregion (IUCN)",
       x = "Mean Fishing Hours (log scale)",
       y = "Mean Conservation Priority (IUCN)") +
  theme_minimal() +
  guides(color = guide_legend(title = "Realm", override.aes = list(size = 4)))

# FUSE - Get top 5 fishing and top 5 priority ecoregions
top_fishing_fuse <- ecoregion_combined_fuse %>%
  arrange(desc(mean_fishing_hours)) %>%
  head(5) %>%
  pull(ECOREGION)

top_priority_fuse <- ecoregion_combined_fuse %>%
  arrange(desc(mean_priority)) %>%
  head(5) %>%
  pull(ECOREGION)

ecoregions_to_label_fuse <- unique(c(top_fishing_fuse, top_priority_fuse))
label_data_fuse <- ecoregion_combined_fuse %>%
  filter(ECOREGION %in% ecoregions_to_label_fuse)

# Create FUSE plot
fuse_ecoregion_plot <- ggplot(ecoregion_combined_fuse,
                             aes(x = mean_fishing_hours, y = mean_priority)) +
  geom_point(aes(size = fishing_cells, color = REALM), alpha = 0.4) +
  geom_label_repel(
    data = label_data_fuse,
    aes(label = ECOREGION),
    size = 3,
    max.overlaps = 10,
    box.padding = 0.5
  ) +
  scale_x_log10() +
  scale_size_continuous(name = "Number of Cells", range = c(1, 8)) +
  scale_color_manual(values = c("red", "blue", "darkgreen", "purple", "orange", "brown", 
                               "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +
  labs(title = "Relationship Between Fishing Pressure and Conservation Priority",
       subtitle = "By Ecoregion (FUSE)",
       x = "Mean Fishing Hours (log scale)",
       y = "Mean Conservation Priority (FUSE)") +
  theme_minimal() +
  guides(color = guide_legend(title = "Realm", override.aes = list(size = 4)))

# EDGE2 - Get top 5 fishing and top 5 priority ecoregions
top_fishing_edge2 <- ecoregion_combined_edge2 %>%
  arrange(desc(mean_fishing_hours)) %>%
  head(5) %>%
  pull(ECOREGION)

top_priority_edge2 <- ecoregion_combined_edge2 %>%
  arrange(desc(mean_priority)) %>%
  head(5) %>%
  pull(ECOREGION)

ecoregions_to_label_edge2 <- unique(c(top_fishing_edge2, top_priority_edge2))
label_data_edge2 <- ecoregion_combined_edge2 %>%
  filter(ECOREGION %in% ecoregions_to_label_edge2)

# Create EDGE2 plot
edge2_ecoregion_plot <- ggplot(ecoregion_combined_edge2,
                              aes(x = mean_fishing_hours, y = mean_priority)) +
  geom_point(aes(size = fishing_cells, color = REALM), alpha = 0.4) +
  geom_label_repel(
    data = label_data_edge2,
    aes(label = ECOREGION),
    size = 3,
    max.overlaps = 10,
    box.padding = 0.5
  ) +
  scale_x_log10() +
  scale_size_continuous(name = "Number of Cells", range = c(1, 8)) +
  scale_color_manual(values = c("red", "blue", "darkgreen", "purple", "orange", "brown", 
                               "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +
  labs(title = "Relationship Between Fishing Pressure and Conservation Priority",
       subtitle = "By Ecoregion (EDGE2)",
       x = "Mean Fishing Hours (log scale)",
       y = "Mean Conservation Priority (EDGE2)") +
  theme_minimal() +
  guides(color = guide_legend(title = "Realm", override.aes = list(size = 4)))

# Create the plots without titles and legends (we'll use a shared legend)

# IUCN Plot
iucn_ecoregion_plot <- ggplot(ecoregion_combined_iucn,
                             aes(x = mean_fishing_hours, y = mean_priority)) +
  geom_point(aes(size = fishing_cells, color = REALM), alpha = 0.4) +
  geom_label_repel(
    data = label_data_iucn,
    aes(label = ECOREGION),
    size = 3,
    max.overlaps = 10,
    box.padding = 0.5
  ) +
  scale_x_log10() +
  scale_size_continuous(name = "Number of Cells", range = c(1, 8)) +
  scale_color_manual(values = c("red", "blue", "darkgreen", "purple", "orange", "brown", 
                               "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +
  labs(x = "Mean Fishing Hours (log scale)",
       y = "Mean Conservation Priority") +
  theme_minimal() +
  theme(
    legend.position = "none",  # Remove legend from individual plots
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA)
  ) +
  guides(
    color = guide_legend(title = "Realm", override.aes = list(size = 4, alpha = 1)),
    size = guide_legend(title = "Number of Cells")
  )

# FUSE Plot
fuse_ecoregion_plot <- ggplot(ecoregion_combined_fuse,
                             aes(x = mean_fishing_hours, y = mean_priority)) +
  geom_point(aes(size = fishing_cells, color = REALM), alpha = 0.4) +
  geom_label_repel(
    data = label_data_fuse,
    aes(label = ECOREGION),
    size = 3,
    max.overlaps = 10,
    box.padding = 0.5
  ) +
  scale_x_log10() +
  scale_size_continuous(name = "Number of Cells", range = c(1, 8)) +
  scale_color_manual(values = c("red", "blue", "darkgreen", "purple", "orange", "brown", 
                               "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +
  labs(x = "Mean Fishing Hours (log scale)",
       y = "Mean Conservation Priority") +
  theme_minimal() +
  theme(
    legend.position = "none",  # Remove legend from individual plots
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA)
  ) +
  guides(
    color = guide_legend(title = "Realm", override.aes = list(size = 4, alpha = 1)),
    size = guide_legend(title = "Number of Cells")
  )

# EDGE2 Plot
edge2_ecoregion_plot <- ggplot(ecoregion_combined_edge2,
                              aes(x = mean_fishing_hours, y = mean_priority)) +
  geom_point(aes(size = fishing_cells, color = REALM), alpha = 0.4) +
  geom_label_repel(
    data = label_data_edge2,
    aes(label = ECOREGION),
    size = 3,
    max.overlaps = 10,
    box.padding = 0.5
  ) +
  scale_x_log10() +
  scale_size_continuous(name = "Number of Cells", range = c(1, 8)) +
  scale_color_manual(values = c("red", "blue", "darkgreen", "purple", "orange", "brown", 
                               "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +
  labs(x = "Mean Fishing Hours (log scale)",
       y = "Mean Conservation Priority") +
  theme_minimal() +
  theme(
    legend.position = "none",  # Remove legend from individual plots
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA)
  ) +
  guides(
    color = guide_legend(title = "Realm", override.aes = list(size = 4, alpha = 1)),
    size = guide_legend(title = "Number of Cells")
  )

# Create one plot with legend for extracting the common legend
legend_plot <- ggplot(ecoregion_combined_iucn,
                     aes(x = mean_fishing_hours, y = mean_priority)) +
  geom_point(aes(size = fishing_cells, color = REALM), alpha = 0.4) +
  scale_size_continuous(name = "Number of Cells", range = c(1, 8)) +
  scale_color_manual(values = c("red", "blue", "darkgreen", "purple", "orange", "brown", 
                               "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +
  guides(
    color = guide_legend(title = "Realm", override.aes = list(size = 4, alpha = 1)),
    size = guide_legend(title = "Number of Cells")
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA)
  )

# Create a blank plot for the legend position
blank_plot <- ggplot() + theme_void()

# Extract the legend as a separate grob
legend_grob <- cowplot::get_legend(legend_plot)

# Use ggpubr::ggarrange to combine plots in a 2x2 grid with legend in position 4
library(ggpubr)
combined_ecoregion_plots <- ggarrange(
  iucn_ecoregion_plot, fuse_ecoregion_plot, 
  edge2_ecoregion_plot, as_ggplot(legend_grob),
  labels = c("(A)", "(B)", "(C)", ""),  # No label for legend position
  ncol = 2, nrow = 2
)

# Ensure white background for the entire combined plot
combined_ecoregion_plots <- combined_ecoregion_plots + 
  theme(plot.background = element_rect(fill = "white", color = NA))

# Print the combined plot
print(combined_ecoregion_plots)

# Save the plot with adjusted dimensions for 2x2 layout and explicit white background
ggsave(here::here("outputs", "ecoregion_fishing_priority_relationship.png"), 
       combined_ecoregion_plots, 
       width = 12, height = 10, dpi = 300, bg = "white")

# Print summary statistics
cat("Summary of ecoregion-level analysis:\n")

Summary of ecoregion-level analysis:

cat("Ecoregions analyzed - IUCN:", nrow(ecoregion_combined_iucn), "\n")

Ecoregions analyzed - IUCN: 229 

cat("Ecoregions analyzed - FUSE:", nrow(ecoregion_combined_fuse), "\n")

Ecoregions analyzed - FUSE: 229 

cat("Ecoregions analyzed - EDGE2:", nrow(ecoregion_combined_edge2), "\n\n")

Ecoregions analyzed - EDGE2: 229 

# Show labeled ecoregions for each index
cat("Ecoregions labeled in IUCN plot:\n")

Ecoregions labeled in IUCN plot:

print(label_data_iucn %>% dplyr::select(ECOREGION, REALM, mean_fishing_hours, mean_priority) %>% arrange(desc(mean_priority)))

# A tibble: 10 × 4
   ECOREGION                              REALM mean_fishing_hours mean_priority
   <chr>                                  <chr>              <dbl>         <dbl>
 1 Torres Strait Northern Great Barrier … Cent…               846.        0.568 
 2 Floridian                              Trop…              2893.        0.522 
 3 Exmouth to Broome                      Cent…              1020.        0.487 
 4 Sunda Shelf/Java Sea                   Cent…              2815.        0.477 
 5 Uruguay-Buenos Aires Shelf             Temp…              3247.        0.439 
 6 Southern China                         Cent…            171703.        0.291 
 7 East China Sea                         Temp…             55657.        0.208 
 8 Gulf of Tonkin                         Cent…             36383.        0.199 
 9 Adriatic Sea                           Temp…             58855.        0.0529
10 Yellow Sea                             Temp…            113047.        0     

cat("\nEcoregions labeled in FUSE plot:\n")

Ecoregions labeled in FUSE plot:

print(label_data_fuse %>% dplyr::select(ECOREGION, REALM, mean_fishing_hours, mean_priority) %>% arrange(desc(mean_priority)))

# A tibble: 9 × 4
  ECOREGION             REALM                   mean_fishing_hours mean_priority
  <chr>                 <chr>                                <dbl>         <dbl>
1 East China Sea        Temperate Northern Pac…             55657.        0.539 
2 Southern Vietnam      Central Indo-Pacific                 5063.        0.443 
3 Ionian Sea            Temperate Northern Atl…             10660.        0.433 
4 Celtic Seas           Temperate Northern Atl…             15475.        0.421 
5 Western Mediterranean Temperate Northern Atl…             14785.        0.419 
6 Southern China        Central Indo-Pacific               171703.        0.310 
7 Yellow Sea            Temperate Northern Pac…            113047.        0.297 
8 Gulf of Tonkin        Central Indo-Pacific                36383.        0.202 
9 Adriatic Sea          Temperate Northern Atl…             58855.        0.0361

cat("\nEcoregions labeled in EDGE2 plot:\n")

Ecoregions labeled in EDGE2 plot:

print(label_data_edge2 %>% dplyr::select(ECOREGION, REALM, mean_fishing_hours, mean_priority) %>% arrange(desc(mean_priority)))

# A tibble: 10 × 4
   ECOREGION                              REALM mean_fishing_hours mean_priority
   <chr>                                  <chr>              <dbl>         <dbl>
 1 Faroe Plateau                          Temp…              6820.         0.729
 2 Torres Strait Northern Great Barrier … Cent…               846.         0.545
 3 Exmouth to Broome                      Cent…              1020.         0.482
 4 Arafura Sea                            Cent…              2336.         0.450
 5 Uruguay-Buenos Aires Shelf             Temp…              3247.         0.443
 6 Adriatic Sea                           Temp…             58855.         0.335
 7 Southern China                         Cent…            171703.         0.295
 8 East China Sea                         Temp…             55657.         0.191
 9 Gulf of Tonkin                         Cent…             36383.         0.153
10 Yellow Sea                             Temp…            113047.         0    

# Optional: Create individual plots if needed
# print(iucn_ecoregion_plot)
# print(fuse_ecoregion_plot) 
# print(edge2_ecoregion_plot)

# The main combined plot is printed above

Maps of SR, FUn and EDGE2

# Sepcies richness ----
# Load required libraries
library(here)
library(terra)
library(ggplot2)
library(viridis)
library(sf)
library(rnaturalearth)
library(smoothr)

# Set the path to your raster files
raster_path <- here("Data", "tif files continental")

# Get list of all TIF files in the directory
tif_files <- list.files(raster_path, pattern = "\\.tif$", full.names = TRUE)

# Check if files were found
if(length(tif_files) == 0) {
  stop("No TIF files found in the specified directory")
}

print(paste("Found", length(tif_files), "TIF files"))

[1] "Found 1000 TIF files"

# Read all rasters into a SpatRaster stack
species_rasters <- rast(tif_files)

# Calculate species richness by summing across all layers
# This assumes presence = 1, absence = 0 or NA
species_richness <- sum(species_rasters, na.rm = TRUE)

# Name the layer
names(species_richness) <- "Species_Richness"

# Create a basic plot
plot(species_richness, 
     main = "Species Richness Map",
     col = viridis(100))

# Alternative: Create a ggplot map with McBryde-Thomas 2 projection
# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Project the species richness raster to McBryde-Thomas 2
species_richness_projected <- project(species_richness, mcbryde_thomas_2)

# Convert to dataframe for ggplot
richness_df <- as.data.frame(species_richness_projected, xy = TRUE)

# Remove NA values for cleaner plotting
richness_df <- richness_df[!is.na(richness_df$Species_Richness), ]

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# Convert richness data to sf points for proper clipping
richness_sf <- st_as_sf(richness_df, coords = c("x", "y"), crs = mcbryde_thomas_2)

# Clip the richness data to the globe boundary
richness_clipped <- st_intersection(richness_sf, globe_border)

# Extract coordinates back for plotting
richness_clipped_coords <- st_coordinates(richness_clipped)
richness_final <- data.frame(
  x = richness_clipped_coords[,1],
  y = richness_clipped_coords[,2],
  Species_Richness = richness_clipped$Species_Richness
)

# Create ggplot map
species_richness_plot <- ggplot() +
  geom_tile(data = richness_final, aes(x = x, y = y, fill = Species_Richness)) +
  geom_sf(data = world_projected, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_fill_viridis_c(
    name = "Species\nRichness",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  coord_sf(crs = mcbryde_thomas_2, expand = FALSE) +
  labs(title = "Global Species Richness Map",
       x = NULL, y = NULL) +
  my_theme

# Display the plot
print(species_richness_plot)

# Richness of threatened species, top 25% FUSE and top 25% EDGE2 ----

# Load required libraries
library(here)
library(terra)
library(ggplot2)
library(viridis)
library(sf)
library(rnaturalearth)
library(smoothr)
library(dplyr)

# Load the conservation metrics dataset
conservation_data <- read.csv(here("Data", "My dataframes", "continental_shark_conservation_metrics_10_harmonised_IUCN_csv.csv"))

# Filter for threatened species (IUCN values of 0.5, 0.75, or 1)
threatened_species <- conservation_data %>%
  filter(IUCN %in% c(0.5, 0.75, 1)) %>%
  pull(Species.name)

# Get top 25% FUSE scoring species by rank (not by threshold)
n_species <- nrow(conservation_data)
top_25_percent_count <- ceiling(n_species * 0.25)

top_fuse_species <- conservation_data %>%
  arrange(desc(FUSE)) %>%
  slice_head(n = top_25_percent_count) %>%
  pull(Species.name)

# Get top 25% EDGE2 scoring species by rank
top_edge2_species <- conservation_data %>%
  arrange(desc(EDGE2)) %>%
  slice_head(n = top_25_percent_count) %>%
  pull(Species.name)

# Get the actual threshold values for reporting
fuse_threshold <- conservation_data %>%
  arrange(desc(FUSE)) %>%
  slice(top_25_percent_count) %>%
  pull(FUSE)

edge2_threshold <- conservation_data %>%
  arrange(desc(EDGE2)) %>%
  slice(top_25_percent_count) %>%
  pull(EDGE2)

print(paste("Total species in dataset:", n_species))

[1] "Total species in dataset: 1000"

print(paste("Top 25% count:", top_25_percent_count))

[1] "Top 25% count: 250"

print(paste("Found", length(threatened_species), "threatened species"))

[1] "Found 366 threatened species"

print(paste("Found", length(top_fuse_species), "top 25% FUSE species (threshold:", round(fuse_threshold, 4), ")"))

[1] "Found 250 top 25% FUSE species (threshold: 0.1 )"

print(paste("Found", length(top_edge2_species), "top 25% EDGE2 species (threshold:", round(edge2_threshold, 4), ")"))

[1] "Found 250 top 25% EDGE2 species (threshold: 0.1 )"

# Set the path to your raster files
raster_path <- here("Data", "tif files continental")

# Get list of all TIF files in the directory
all_tif_files <- list.files(raster_path, pattern = "\\.tif$", full.names = TRUE)

# Extract species names from file names (assuming filename format includes species name)
# You may need to adjust this pattern based on your actual file naming convention
file_species_names <- basename(all_tif_files)
file_species_names <- gsub("\\.tif$", "", file_species_names)

# Function to find matching raster files for a species list
find_matching_files <- function(species_list, all_files) {
  matching_files <- c()
  for(species in species_list) {
    # Create pattern to match species name in filename (handles spaces and underscores)
    species_pattern <- gsub(" ", "[_ ]", species)
    matches <- all_files[grepl(species_pattern, basename(all_files), ignore.case = TRUE)]
    matching_files <- c(matching_files, matches)
  }
  return(unique(matching_files))
}

# Find matching files for each species group
threatened_files <- find_matching_files(threatened_species, all_tif_files)
fuse_files <- find_matching_files(top_fuse_species, all_tif_files)
edge2_files <- find_matching_files(top_edge2_species, all_tif_files)

print(paste("Found", length(threatened_files), "raster files for threatened species"))

[1] "Found 366 raster files for threatened species"

print(paste("Found", length(fuse_files), "raster files for top FUSE species"))

[1] "Found 250 raster files for top FUSE species"

print(paste("Found", length(edge2_files), "raster files for top EDGE2 species"))

[1] "Found 250 raster files for top EDGE2 species"

# Check if files were found
if(length(threatened_files) == 0) {
  stop("No matching TIF files found for threatened species. Please check file naming convention.")
}

# Read rasters and calculate richness for each group
# Threatened species
threatened_rasters <- rast(threatened_files)
threatened_richness <- sum(threatened_rasters, na.rm = TRUE)
names(threatened_richness) <- "Threatened_Species_Richness"

# Top FUSE species
fuse_rasters <- rast(fuse_files)
fuse_richness <- sum(fuse_rasters, na.rm = TRUE)
names(fuse_richness) <- "Top_FUSE_Richness"

# Top EDGE2 species
edge2_rasters <- rast(edge2_files)
edge2_richness <- sum(edge2_rasters, na.rm = TRUE)
names(edge2_richness) <- "Top_EDGE2_Richness"

# Create basic plots for all three metrics
plot(threatened_richness, 
     main = "Threatened Species Richness Map",
     col = viridis(100))

plot(fuse_richness, 
     main = "Top 25% FUSE Species Richness Map",
     col = viridis(100))

plot(edge2_richness, 
     main = "Top 25% EDGE2 Species Richness Map",
     col = viridis(100))

# Function to create projected richness map
create_richness_map <- function(richness_raster, title, subtitle, legend_name) {
  # Project the richness raster to McBryde-Thomas 2
  richness_projected <- project(richness_raster, mcbryde_thomas_2)
  
  # Convert to dataframe for ggplot
  richness_df <- as.data.frame(richness_projected, xy = TRUE)
  colnames(richness_df)[3] <- "richness_value"
  
  # Remove NA values for cleaner plotting
  richness_df <- richness_df[!is.na(richness_df$richness_value), ]
  
  # Convert richness data to sf points for proper clipping
  richness_sf <- st_as_sf(richness_df, coords = c("x", "y"), crs = mcbryde_thomas_2)
  
  # Clip the richness data to the globe boundary
  richness_clipped <- st_intersection(richness_sf, globe_border)
  
  # Extract coordinates back for plotting
  richness_clipped_coords <- st_coordinates(richness_clipped)
  richness_final <- data.frame(
    x = richness_clipped_coords[,1],
    y = richness_clipped_coords[,2],
    richness_value = richness_clipped$richness_value
  )
  
  # Create ggplot map
  richness_plot <- ggplot() +
    geom_tile(data = richness_final, aes(x = x, y = y, fill = richness_value)) +
    geom_sf(data = world_projected, fill = "lightgrey", color = "lightgrey", size = 0.1) +
    geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
    scale_fill_viridis_c(
      name = legend_name,
      guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                           title.position = "top", title.hjust = 0.5)
    ) +
    coord_sf(crs = mcbryde_thomas_2, expand = FALSE) +
    labs(title = title,
         subtitle = subtitle,
         x = NULL, y = NULL) +
    my_theme
  
  return(richness_plot)
}

# Create maps for all three metrics
threatened_richness_plot <- create_richness_map(
  threatened_richness, 
  "Global Threatened Species Richness Map",
  "Species with IUCN threat levels 0.5, 0.75, or 1.0",
  "Threatened\nSpecies Richness"
)

fuse_richness_plot <- create_richness_map(
  fuse_richness,
  "Global Top 25% FUSE Species Richness Map", 
  paste("Species in top quartile of FUSE scores (≥", round(fuse_threshold, 4), ")"),
  "Top FUSE\nSpecies Richness"
)

edge2_richness_plot <- create_richness_map(
  edge2_richness,
  "Global Top 25% EDGE2 Species Richness Map",
  paste("Species in top quartile of EDGE2 scores (≥", round(edge2_threshold, 4), ")"),
  "Top EDGE2\nSpecies Richness"
)

# Display all plots
print(threatened_richness_plot)

print(fuse_richness_plot)

print(edge2_richness_plot)

Differences between priority and normal richness maps

create_comparison_maps <- function(captain2_data, richness_raster, priority_threshold = 0.85, 
                                   comparison_name = "Comparison") {
  
  # Count grid cells with CAPTAIN2 priority > threshold
  high_priority_cells <- sum(captain2_data$Priority > priority_threshold, na.rm = TRUE)
  
  print(paste("CAPTAIN2", comparison_name, "cells with priority >", priority_threshold, ":", high_priority_cells))
  
  # Get the same number of top richness cells
  richness_values <- values(richness_raster, na.rm = TRUE)
  richness_threshold <- quantile(richness_values, 1 - (high_priority_cells / length(richness_values)), na.rm = TRUE)
  
  print(paste("Richness threshold for top", high_priority_cells, "cells:", round(richness_threshold, 2)))
  
  # Create binary rasters for comparison
  # CAPTAIN2 binary (1 = high priority, 0 = other)
  captain2_binary <- captain2_data
  captain2_binary$Binary_Priority <- ifelse(captain2_data$Priority > priority_threshold, 1, 0)
  
  # Richness binary (1 = top cells, 0 = other)  
  richness_binary <- richness_raster
  values(richness_binary) <- ifelse(values(richness_raster) >= richness_threshold, 1, 0)
  
  # Convert CAPTAIN2 data to sf if it's not already
  if(!inherits(captain2_binary, "sf")) {
    # Assuming captain2_data has Longitude and Latitude columns
    captain2_binary_sf <- st_as_sf(captain2_binary, 
                                   coords = c("Longitude", "Latitude"), 
                                   crs = 4326)
  } else {
    captain2_binary_sf <- captain2_binary
  }
  
  # Transform to match the richness raster CRS
  captain2_binary_sf <- st_transform(captain2_binary_sf, crs = crs(richness_binary))
  
  # Convert to raster using terra::rasterize with explicit method
  captain2_raster <- rasterize(captain2_binary_sf, richness_binary, 
                               field = "Binary_Priority", 
                               fun = "max",  # Use max to avoid averaging
                               background = 0)  # Set background to 0
  
  # Ensure we only have 0s and 1s
  captain2_raster[captain2_raster > 0] <- 1
  
  # Define the McBryde-Thomas 2 projection
  mcbryde_thomas_2 <- "+proj=mbt_s"
  
  # Calculate difference and debug each step
  difference_raster <- captain2_raster - richness_binary
  print(paste("Difference raster total cells:", ncell(difference_raster)))
  print(paste("Difference raster non-NA cells:", sum(!is.na(values(difference_raster)))))
  
  # Create a mask for cells where at least one method has priority (value = 1)
  priority_mask <- (captain2_raster == 1) | (richness_binary == 1)
  print(paste("Priority mask TRUE cells:", sum(values(priority_mask), na.rm = TRUE)))
  
  # Apply the mask to keep only priority cells
  difference_raster[!priority_mask] <- NA
  print(paste("After masking, non-NA cells:", sum(!is.na(values(difference_raster)))))
  
  # Project difference raster
  difference_projected <- project(difference_raster, mcbryde_thomas_2)
  print(paste("After projection, non-NA cells:", sum(!is.na(values(difference_projected)))))
  
  # Convert to dataframe and clip to globe
  diff_df <- as.data.frame(difference_projected, xy = TRUE)
  colnames(diff_df)[3] <- "difference"
  diff_df <- diff_df[!is.na(diff_df$difference), ]
  print(paste("After df conversion, rows:", nrow(diff_df)))
  
  # Convert to sf and clip
  diff_sf <- st_as_sf(diff_df, coords = c("x", "y"), crs = mcbryde_thomas_2)
  print(paste("SF object rows:", nrow(diff_sf)))
  
  # Create world map and globe border for this function
  world <- ne_countries(scale = "medium", returnclass = "sf")
  world_projected <- st_transform(world, crs = mcbryde_thomas_2)
  
  # Create the globe bounding box and border
  globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                      c(180, 90), c(180, -90), c(-180, -90))
  globe_border <- st_polygon(list(globe_bbox)) %>%
    st_sfc(crs = 4326) %>%
    st_sf(data.frame(rgn = 'globe', geom = .)) %>%
    smoothr::densify(max_distance = 0.5) %>%
    st_transform(crs = mcbryde_thomas_2)
  
  # Create base theme
  my_theme <- theme_minimal() +
    theme(
      legend.position = "bottom",
      legend.direction = "horizontal",
      legend.box = "vertical",
      legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
      legend.title = element_text(margin = margin(b = 10)),
      panel.background = element_rect(fill = "white", color = NA),
      plot.background = element_rect(fill = "white", color = NA),
      panel.grid = element_blank()
    )
  
  diff_clipped <- st_intersection(diff_sf, globe_border)
  print(paste("After globe clipping, rows:", nrow(diff_clipped)))
  
  # Extract coordinates
  diff_coords <- st_coordinates(diff_clipped)
  diff_final <- data.frame(
    x = diff_coords[,1],
    y = diff_coords[,2],
    difference = diff_clipped$difference
  )
  print(paste("Final data rows:", nrow(diff_final)))
  
  # Create difference map
  diff_plot <- ggplot() +
    geom_tile(data = diff_final, aes(x = x, y = y, fill = factor(difference))) +
    geom_sf(data = world_projected, fill = "lightgrey", color = "lightgrey", size = 0.1) +
    geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
    scale_fill_manual(
      values = c("-1" = "red", "0" = "black", "1" = "blue"),
      labels = c("-1" = "Only Richness", "0" = "Agreement", "1" = "Only CAPTAIN2"),
      name = "Priority\nDifference",
      guide = guide_legend(title.position = "top", title.hjust = 0.5)
    ) +
    coord_sf(crs = mcbryde_thomas_2, expand = FALSE) +
    labs(title = paste("Priority Comparison:", comparison_name),
         subtitle = paste("CAPTAIN2 vs Simple Richness (", high_priority_cells, "top cells each)"),
         x = NULL, y = NULL) +
    my_theme +
    theme(legend.position = "bottom")
  
  # Calculate summary statistics
  total_cells <- nrow(diff_final)
  agreement_cells <- sum(diff_final$difference == 0, na.rm = TRUE)
  captain2_only <- sum(diff_final$difference == 1, na.rm = TRUE)
  richness_only <- sum(diff_final$difference == -1, na.rm = TRUE)
  
  agreement_percent <- (agreement_cells / total_cells) * 100
  captain2_only_percent <- (captain2_only / total_cells) * 100
  richness_only_percent <- (richness_only / total_cells) * 100
  
  # Print summary statistics with debugging info
  cat("\n========================================\n")
  cat("SUMMARY STATISTICS -", comparison_name, "\n")
  cat("========================================\n")
  cat("CAPTAIN2 priority cells:", high_priority_cells, "\n")
  cat("Richness priority cells:", high_priority_cells, "(same number)\n")
  cat("CAPTAIN2 raster 1s:", sum(values(captain2_raster) == 1, na.rm = TRUE), "\n")
  cat("Richness raster 1s:", sum(values(richness_binary) == 1, na.rm = TRUE), "\n")
  cat("Total priority cells compared:", total_cells, "\n")
  cat("----------------------------------------\n")
  cat("Agreement (both prioritize same cell):", agreement_cells, "(", round(agreement_percent, 2), "%)\n")
  cat("Only CAPTAIN2 prioritizes:", captain2_only, "(", round(captain2_only_percent, 2), "%)\n")
  cat("Only Richness prioritizes:", richness_only, "(", round(richness_only_percent, 2), "%)\n")
  cat("----------------------------------------\n")
  cat("Total accounted for:", round(agreement_percent + captain2_only_percent + richness_only_percent, 2), "%\n")
  cat("========================================\n\n")
  
  return(list(plot = diff_plot, 
              captain2_cells = high_priority_cells,
              richness_threshold = richness_threshold,
              total_cells = total_cells,
              agreement_cells = agreement_cells,
              agreement_percent = agreement_percent,
              captain2_only_percent = captain2_only_percent,
              richness_only_percent = richness_only_percent))
}

# Threatened species comparison
threatened_comparison <- create_comparison_maps(
  CAPTAIN2_IUCN_data_nonzero, 
  threatened_richness, 
  priority_threshold = 0.85,
  comparison_name = "Threatened Species"
)

[1] "CAPTAIN2 Threatened Species cells with priority > 0.85 : 4875"
[1] "Richness threshold for top 4875 cells: 25"
[1] "Difference raster total cells: 232560"
[1] "Difference raster non-NA cells: 47224"
[1] "Priority mask TRUE cells: 6203"
[1] "After masking, non-NA cells: 6203"
[1] "After projection, non-NA cells: 38272"
[1] "After df conversion, rows: 38272"
[1] "SF object rows: 38272"

[1] "After globe clipping, rows: 37964"
[1] "Final data rows: 37964"

========================================
SUMMARY STATISTICS - Threatened Species 
========================================
CAPTAIN2 priority cells: 4875 
Richness priority cells: 4875 (same number)
CAPTAIN2 raster 1s: 4875 
Richness raster 1s: 5180 
Total priority cells compared: 37964 
----------------------------------------
Agreement (both prioritize same cell): 20049 ( 52.81 %)
Only CAPTAIN2 prioritizes: 4872 ( 12.83 %)
Only Richness prioritizes: 5838 ( 15.38 %)
----------------------------------------
Total accounted for: 81.02 %
========================================

print(threatened_comparison$plot)

# FUSE comparison  
fuse_comparison <- create_comparison_maps(
  CAPTAIN2_FUSE_data_nonzero,
  fuse_richness,
  priority_threshold = 0.85, 
  comparison_name = "FUSE"
)

[1] "CAPTAIN2 FUSE cells with priority > 0.85 : 4601"
[1] "Richness threshold for top 4601 cells: 17"
[1] "Difference raster total cells: 232560"
[1] "Difference raster non-NA cells: 49658"
[1] "Priority mask TRUE cells: 6964"
[1] "After masking, non-NA cells: 6964"
[1] "After projection, non-NA cells: 41088"
[1] "After df conversion, rows: 41088"
[1] "SF object rows: 41088"

[1] "After globe clipping, rows: 40782"
[1] "Final data rows: 40782"

========================================
SUMMARY STATISTICS - FUSE 
========================================
CAPTAIN2 priority cells: 4601 
Richness priority cells: 4601 (same number)
CAPTAIN2 raster 1s: 4601 
Richness raster 1s: 4790 
Total priority cells compared: 40782 
----------------------------------------
Agreement (both prioritize same cell): 10036 ( 24.61 %)
Only CAPTAIN2 prioritizes: 9271 ( 22.73 %)
Only Richness prioritizes: 11962 ( 29.33 %)
----------------------------------------
Total accounted for: 76.67 %
========================================

print(fuse_comparison$plot)

# EDGE2 comparison
edge2_comparison <- create_comparison_maps(
  CAPTAIN2_EDGE2_data_nonzero,
  edge2_richness,
  priority_threshold = 0.85,
  comparison_name = "EDGE2" 
)

[1] "CAPTAIN2 EDGE2 cells with priority > 0.85 : 5341"
[1] "Richness threshold for top 5341 cells: 14"
[1] "Difference raster total cells: 232560"
[1] "Difference raster non-NA cells: 48863"
[1] "Priority mask TRUE cells: 7515"
[1] "After masking, non-NA cells: 7515"
[1] "After projection, non-NA cells: 45346"
[1] "After df conversion, rows: 45346"
[1] "SF object rows: 45346"

[1] "After globe clipping, rows: 43802"
[1] "Final data rows: 43802"

========================================
SUMMARY STATISTICS - EDGE2 
========================================
CAPTAIN2 priority cells: 5341 
Richness priority cells: 5341 (same number)
CAPTAIN2 raster 1s: 5341 
Richness raster 1s: 5436 
Total priority cells compared: 43802 
----------------------------------------
Agreement (both prioritize same cell): 15169 ( 34.63 %)
Only CAPTAIN2 prioritizes: 8965 ( 20.47 %)
Only Richness prioritizes: 11058 ( 25.25 %)
----------------------------------------
Total accounted for: 80.34 %
========================================

print(edge2_comparison$plot)